

%% load estimation results
load('parameter_estimates_direction_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
KZ=size(Z,2)/5;
nBeta = length(reb{1});
nPhi1 = length(mpow{1});
nPhi2 = length(mpow{2});
nGamma1 = length(reb{2});
nGamma2 = length(reb{3});

gamma_US1=thetaStar(nBeta+4+nPhi1+4+nPhi2+1,1);
gamma_UK1=thetaStar(nBeta+4+nPhi1+4+nPhi2+2,1);
gamma_FRN1=thetaStar(nBeta+4+nPhi1+4+nPhi2+3,1);
gamma_RUS1=thetaStar(nBeta+4+nPhi1+4+nPhi2+4,1);
gamma_CHN1=thetaStar(nBeta+4+nPhi1+4+nPhi2+5,1);
gamma01=thetaStar(nBeta+4+nPhi1+4+nPhi2+5+(1:nGamma1),1);
gamma_US2=thetaStar(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+1,1);
gamma_UK2=thetaStar(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+2,1);
gamma_FRN2=thetaStar(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+3,1);
gamma_RUS2=thetaStar(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+4,1);
gamma_CHN2=thetaStar(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5,1);
gamma02=thetaStar(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+(1:nGamma2),1);

ZUS = Z(:,1:KZ);
ZUK = Z(:,KZ+(1:KZ)); 
ZFRN = Z(:,2*KZ+(1:KZ));
ZRUS = Z(:,3*KZ+(1:KZ));
ZCHN = Z(:,4*KZ+(1:KZ));

ZrebelUS1 = ZUS(:,reb{2}); ZrebelUS2 = ZUS(:,reb{3});
ZrebelUK1 = ZUK(:,reb{2}); ZrebelUK2 = ZUK(:,reb{3});
ZrebelFRN1 = ZFRN(:,reb{2}); ZrebelFRN2 = ZFRN(:,reb{3});
ZrebelRUS1 = ZRUS(:,reb{2}); ZrebelRUS2 = ZRUS(:,reb{3});
ZrebelCHN1 = ZCHN(:,reb{2}); ZrebelCHN2 = ZCHN(:,reb{3});


%% compute effect of gov-sided intervention
M1 = [mean(gamma_US1+ZrebelUS1*gamma01),...
    mean(gamma_UK1+ZrebelUK1*gamma01),...
    mean(gamma_FRN1+ZrebelFRN1*gamma01),...
    mean(gamma_RUS1+ZrebelRUS1*gamma01),...
    mean(gamma_CHN1+ZrebelCHN1*gamma01)];

S1 = [std(gamma_US1+ZrebelUS1*gamma01),...
    std(gamma_UK1+ZrebelUK1*gamma01),...
    std(gamma_FRN1+ZrebelFRN1*gamma01),...
    std(gamma_RUS1+ZrebelRUS1*gamma01),...
    std(gamma_CHN1+ZrebelCHN1*gamma01)];


%% compute effect of reb-sided intervention
M2 = [mean(gamma_US2+ZrebelUS2*gamma02),...
    mean(gamma_UK2+ZrebelUK2*gamma02),...
    mean(gamma_FRN2+ZrebelFRN2*gamma02),...
    mean(gamma_RUS2+ZrebelRUS2*gamma02),...
    mean(gamma_CHN2+ZrebelCHN2*gamma02)];

S2 = [std(gamma_US2+ZrebelUS2*gamma02),...
    std(gamma_UK2+ZrebelUK2*gamma02),...
    std(gamma_FRN2+ZrebelFRN2*gamma02),...
    std(gamma_RUS2+ZrebelRUS2*gamma02),...
    std(gamma_CHN2+ZrebelCHN2*gamma02)];


%%
use=[M1; S1; M2; S2]';
use=round(use, 3);
disp(' ')
disp('for Table F2:')
disp(use)



